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ABSTRACT 

We investigate stationary slim accretion disks around Kerr black holes. We construct a 
new numerical method based on the relaxation technique. We systematically cover the whole 
parameter space relevant to stellar mass X-ray binaries. We also notice some non-monotonic 
features in the disk structure, overlooked in previous studies. 
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1. INTRODUCTION 

Obse rvations of microquasars in the rmal states 
(see e.g. McClintock fe Remi llard 20 031) , as well as 
obser vations of some AGN (e.g. lCollin fe Kawaguchi 
2004), point to the existence of black hole (BH) 
accretion flows that are quasi-steady for relatively 
long times, optically thick, and geomerically thin. 
They exist in a rather wide range of accretion 
rates. Many authors study such flows using the 
celebrated Shakura-Sunyaev model. The model 
is easy in applications because of its remarkable 
advantages, in particular: [1] Within the disk, 
matter rotates on almost exact circular Keplerian 
orbits. [2] The inner edge of the disk, i.e. the 
place where the accretion flow changes from rota- 
tion to plunging, always locates at r„ = ISCO. 
[3] The local flux of radiation F(r) emitted from 
the disk surface is given by a universal analytic 
formula, derived directly from conservation of 
mass, energy, and angular momentum, and in- 
dependent of dissipation. The general relativistic 
version of it, F(r) = F nt(t), was derived by 
Novikov fc Thornd l|l973h . 

However, for high accretion rates, e.g. M > 
O.lMsdd, these properties do not hold, and the 
Shakura-Sunyaev model should be changed in or- 
der to properly describe the "slim disk" effects 



(jAbramowicz et alJ [l988h : [1*] Within most of 



the disk the rotation is slightly sub-Keplerian, 
and in the innermost part of the disk the rota- 
tion is super-Keplerian. [2*] With the increas- 
ing accretion rate, the inner edge of the disk 



r in goes closer to the BH than the ISCO (e.g. 
Arj n /nsco ~ 10% for the Eddington mass accre- 
tion rate). [3*] The local flux of radiation emitted 
from the disk surface is not given by the universal 
analytic Novikov-Thorne formula, but should be 
calculated to include advection of heat captured 
within the accreted matter. Slim disks are far less 
radiatively efficient than the Shakura-Sunyaev so- 
lutions. In this respect they resemble advection 
dominated accretion flows (fo r a recent review see 
Narayan fe McClintockl l2008t ) and together with 



them form a class or radiatively inefficient flows 
(RIFs). 

Most of the RIFs models in Kerr geometn 
constructed to 



tuts models m Kerr geometry 
dat e dAbramowicz et al. 1996; 



Gammie fc Pophaml [1998) described advection 
dominated flows rather than slim disks. The 
known slim disk solutions do not provide a suf- 
ficient coverage of the parameter space needed 
for accurate fitting of the calculated continuum 
spectra to these o bserved. The previous fits (e.g. 



IShafee et al. 20061 ) have been done with the as- 
sumption that ri n = ISCO and F(r) — FNT(r)- 
This is clearly ina dequate for L > O .SLEdd, as 
Figures 7 and 8 in Shafee et al. ( 2006h show. In 



this paper, we have calculated a large set of slim 
disk models that cover all the relevant parameter 
space. Hopefully, they will improve the accuracy 
of the spectral fits as well as other astrophysical 
applications of black hole accretion disk models. 

The paper is organized as following: At the be- 
ginning (|2]) we present equations governing rela- 
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tivistic accretion disk for optically thick case. In 
$3]we describe in details the numerical methods we 
use to solve the problem. The solutions for vari- 
ous disk parameters are presented in §31 In §4.21 
we put attention to non-monotonic features in ra- 
dial velocity profiles which appear for some par- 
ticular mass accretion rate range. In the follow- 
ing paragraph (jj5|) we discuss properties of slim 
disk solutions in different regimes of the system 
parameters. Finally, in §7] we summarize our work 
and discuss astrophysical application of the results 
presented here. 

2. RELATIVISTIC SLIM DISKS 

In this section we present slim disk equations 
derived in the Kerr spacetime metric. We fol- 
low the authors who have deriv ed thes e equa- 
tions previously, starting from Lasota (.1994) 
and taking into ac c ount f urthe r improvements by 
Abr amowicz et al. ( 19961 ) and Abramowicz et al. 
( 19971 ). Similar models were constructed and 



solved bv lGammie &: Popham (1998) and B eloborodov 
(1998). We assume: the Kerr metric, axis sym- 
metry (dcf, = 0), stationarity (dt — 0), zero torque 
close to the horizon, large optical depth, no self- 
irradiation, no angular momentum taken away by 
radiation and we neglect the magnetic pressure, 
(i) The mass conservation: 



M 



-2tt£A 



1/2 _ 



V 



(1) 



where £ = f\ pdz is disk surface density and V, 
defined by the relation u r = VA 1 / 2 /(rVT~^V 2 ), 
is the gas radial velocity as measured by an ob- 
server at fixed r who corotates with the fluid. A 
and A are the standard Kerr metric coefficients, 
(ii) The radial momentum conservation: 
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and f2 = u^/it* is the angular velocity with re- 
spect to the stationary observer, f2 = SI — u> is 
the angular velocity with respect to the inertial 
observer, n± = ±M 1 / 2 /(r 3 / 2 ± aM 1 ' 2 ) are the 



angular frequencies of the corotating and counter- 
rotating Keplerian orbits and R = Aj(r 2 A 1 ! 2 ) is 
the radius of gyration. 

(iii) The angular momentum conservation 
(jLasotalll994l ): 



M s A 1 ' 2 A 1 ' 2 1 

— {L-L in ) = ±aP 4 

2ir r 

where L — u<p, £j„ is the angular momentum at 
the disk inner edge, 7 is the Lorentz factor and 
P = 2Hp can be considered vertically integrated 
pressure. Herein, unless stated otherwise, we as- 
sume a = 0.1. 



(iy) The vertical equilibrium (jAbramowicz et al 



1997): 
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with e = ut being the conserved energy associated 
with the Killing vector. 

(v) The energy conservation: 
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where the disk central temperature T has been 
introduced. 

(vi) The regularity condition: 
By a series of algebraic manipulations of the above 
equations we get: 
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with Af and V given by: 
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where : B, C and 1^ are defined in lAbramowicz et alJ 
(1996), P is gas to total pressure ratio and F adv is 
given by Eq. [6l To obtain a physical solution Af 
and V must vanish at the same radius called the 
sonic point. 
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3. NUMERICAL METHODS 



To obtain slim disk solutions one has to solve 
a two dimensional system of ordinary differential 
equations together with the following regularity 
conditions at the sonic radius rs: 







(10) 



and outern boundary conditions given at some 
large radius r out . The location of the sonic point 
is not known a priori. A mathematical problem 
defined in such a way can be classified as a two 
point boundary value problem w hich we soly e ap- 
plying the relaxation technique ijPressl 12002) . To 
start the relaxation process one has to provide a 
trial solution. The convergence strictly depends 
on the quality of such initial guess. In this work 
we apply the following method of searching for the 
proper initial conditions for relaxation. 

The angular momentum at the inner edge of 
the disk, £j n , is the eigen value of the problem 
and must be chosen properly to satisfy the regu- 
larity conditions given by Eq. [lOl The value of Li n 
determines the shape of a solution. It turns out 
that the topology of slim disk solutions with d n 
higher and C m lower than the proper value Ci n ,o 
is different. The former branch of solutions termi- 
nates as soon as the denominator V vanishes (the 
regularity condition is not satisfied). The latter 
does not cross the sonic point at all (V is always 
positive). The self-consistent solution is expected 
to be the common limit of these two branches. 
A trial solution can be achieved in a few itera- 
tion steps[j] We start integration at r > lOOOrs 
assuming the lNovikov fc Thornd l| 1973ft boundary 
conditions and we use the implicit Runge-Kutta 
method of the 4th order. 



lr To save computational time when looking for the trial so- 
lution we estimate ^ using the d iffusive form of viscosity 
(Eq. 35 infAbramowi cz et alll996h instead of calculating it 
numerically - such an approach is precise enough to ensure 
convergence. This fact explains why the trial and relaxed 
models in Fig. [T]do not coincide. 
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Fig. 1. — Temperature profiles in the vicinity of 
the sonic point for a few iterations leading to the 
trial model used as initial condition in the relax- 
ation procedure (for details see $3]). According to 
the topology solutions with too high value of Ci n 
terminate before reaching the proper sonic point 
while solutions with too low £j„ go through the 
sonic point radius but follow an improper branch. 
The relaxed solution and the location of the sonic 
point are also presented. Models calculated for a 
non-spinning 9.4 M & BH with fa = O.lifiEdd- 

Once the trial solution is found one can start 
relaxation process between r out and the estimated 
position of the sonic radius rs. The standard ap- 
proach has to be modified due to the fact that we 
expect singularity at the inner boundary. Thus, 
we treat the problem as a free boundary prob- 
lem and introduce one more var iable descri bing 
the position of the critical point I Press! [2002ft . In 
this work we usually use 100 mesh points spaced 
logarithmically in radius. 

The radial derivatives d In Q /d In r (Eq. [9]) and 
dfl/dr (Eq. [6]) are evaluated numerically basing 
on the Q and fl profiles in the previous iteration 
step. A relaxed solution is obtained in a few itera- 
tion steps and is then used as an initial condition 
for relaxation when looking for the solution of a 
problem with one parameter (e.g. mass accretion 
rate fn or BH spin a*) slightly changed. For new 
system parameters we look for the outern bound- 
ary conditions at r out by integrating the equations 
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from l.lr out (assuming Keplerian disk there) until 
we reach r out . In such a way a full spectrum of 
system parameters can be achieved effectively. 

Once a solution outside the sonic point is found 
we numerically estimate the radial derivatives of 
V and T at the sonic point using values given at 
r > rs- Taking them into account we make a small 
step from the innermost mesh point (which corre- 
sponds to the sonic point) inward. Then we start 
integrating using standard Runge-Kutta method 
until we get close enough to the horizon. 

4. RESULTS 
4.1. Flux profiles 

For very low mass accretion rates most of en- 
ergy generated at each radius is immediately emit- 
ted away and disk is expected to be very thin. 
Therefore, solutions in that re gime are in general 
consistent with solutions of Novikov fc Thornd 
The slim disk solutions deviate from thin 
disk models when advection becomes important. 
According to Fig. [2] in case of a non-rotating BH 
an accretion disk is no longer radiatively efficient 
for mass accretion rates exceeding 0.1 to Edd where 
riiEdd is the critical mass accretion rate defined as: 



rriEdd 



64ttGM 18 M 
= 2.23 x 10 18 - 



M(7 



g-s 



1 (11) 



corresponding to (in case of a non-rotating BH) 
a disk with the Eddington luminosity (see i )4.4p . 
The rate of heat advected increases with mass ac- 
cretion rate. One can expect that for very high 
mass accretio n rates advection would dominate 
heat transfer ('Abram owicz et al. 1983 E It is im- 



portant to note that inside some particular radius 
(e.g. 8M for to = O.GrfiEdd) heat is no longer 
accumulated in the accreted matter but is gradu- 
ally radiated away amplifying the emission com- 
ing out of viscous processes. This fact has strong 
influence on the radial profiles of the flux emit- 
ted from both sides of an accretion disks (Fig. 
[3]). For mass accretion rates implying significant 
amount of advection the maximum of emission 



is shifted inward: from ~ 10M expected for ra- 
diatively efficient disks around non-spinning BH 
down to 7.5M for to = 0.9m em and even further 
for higher accretion rates. The effect is even more 
distinct for rotating BHs (middle and bottom pan- 
els of Fig. 02 please note the common shift due to 
decreasing radius of the innermost stable circular 
orbit). 




Fig. 2. — Ratio of energy advected to generated 
at each radius for different mass accretion rates 
for a* — 0. Positive values denote region where 
heat generated by viscous processes is stored in the 
accreted matter. Negative values mark radii where 
the cumulated heat is radiated away summing up 
to the flux generated by viscous processes. 



2 There is another class of accretion disks which are advection 
dominated and opticaly thin named Advection Dominated 



Accretion Flow s investigated by a number of grou 



Ichimarul dl977T). iRees et all dl982T).lNaravan fc Yi , 
Naravan fc YT fl995l and I AbramowicT*etaT ||1995| 1 



i ds e.g . 
(1994), 
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Fig. 3. — Flux profiles for different mass accretion 
rates and BH spins. Each subplot contains six 
solid lines for the following mass accretion rates: 
0.01 (the thickest line), 0.1, 0.3, 0.6 and 0.9M E dd 
(the thinnest line). The upper panel is for a non 
rotating BH (a* — 0), the middle one for a* = 0.6 
while the bottom one for a highly spinning BH 
(a* = 0.98). BH mass is 9.4 Mq. 

4.2. Disk structure 

Profiles of a few parameters describing an ac- 
cretion disk for a number of mass accretion rates 
in case of a non-spinning BH are presented in Fig. 
HI The profiles of the radial velocity as measured 
in the fluid corotating frame V are presented in 
the upper left panel. The ratio of the radiation to 
gas components of the total vertically integrated 
pressure P are drawn in the upper right corner. 
Disk thickness and central temperatures are pre- 
sented in the lower panels. 

The radial velocity for large radii coincide with 
values given by Novikov & Thorne solutions which 
are assumed as the outer boundary conditions and 
approaches the speed of light when getting close to 
the horizon. At a given radius the radial velocity 



increases with mass accretion rate. For a specific 
range of accretion rates the radial velocity profiles 
are not everywhere monotonic i.e. the velocity is 
not steadily increasing when gas is getting closer 
to the horizon. These features, existence of which 
is out of common knowledge, correspond to regions 
of increased surface density and will be discussed 
in details in the following paragraph. 




Fig. 4. — Profiles of radial velocity (upper-left 
panel), radiation to gas pressure ratio (upper- 
right), disk thickness over radius (H/r) ratio 
(bottom-left) and central temperature (bottom- 
right) for a nonspinning BH of mass 9.4 Mq. The 
solid lines are for mass accretion rates between 
0.01 and 0.9mEdd- 



Shakura &; Sunvaev 1 1973t ) pointed out that an 



accretion disk can be divided into three distinc- 
tive regimes: (i) the outer region, gas pressure 
and free-free scatterings dominated; (ii) the middle 
region, gas pressure but electron scattering dom- 
inated; (iii) the inner region, radiation pressure 
and electron scattering dominated. According to 
their inner region formulae the surface density is 
expected to rise infinitely when approaching the 
inner edge of the disk (accreted gas is slowed down 
by increasing radiation pressure). Obviously, this 
behaviour (which corresponds to a decrease of the 
fluid radial velocity) is suppressed as the Shakura 
& Sunyaev's inner region does not infact extend 
down to the inner edge of the disk. The radiation 
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to gas pr essure ratio in their soluti ons is propor- 
tional to 1 Shakura fc Sunvaev 1973h : 



The transition between inner and middle regions 
takes place where gas and radiation components 
of pressure are equal. It is trivial to show that 
the inner region must be followed by a transition 
to the middle one before reaching the disk inner 
boundary. Therefore, close to 6M the growth of 
the surface density is suppressed by a decrease ex- 
pected in the gas pressure dominated (middle) re- 
gion of an accretion disk. The resulting humps 
in the surface density as well as the correspond- 
ing non-monotonic parts of the radial velocity pro- 
files are presented in details in Fig. \E[ According 
to the upper-right panel of Fig. [4] for the lowest 
mass accretion rates (m < O.OliiiEdd) radiation 
pressure never dominates in the disk and there- 
fore there is no hump in surface density profile at 
all. For moderate mass accretion rates radiation 
pressure exceeds gas pressure in wide range of radii 
and resulting humps in surface density profiles are 
clearly visible. When accretion rate increases the 
advection of heat becomes more and more impor- 
tant (see i )4.ip and the Shakura & Sunyaev for- 
malism cannot be longer applied. Our solutions 
show that for mass accretion rates higher than 
0.9riiEdd (for a* = 0) the surface density and ra- 
dial velocity profiles become again monotonic for 
all radii despite the fact that inner parts of the disk 
are radiation pressure dominated. It is important 
to understand that such humps in surface density 
profiles cannot be considered shock like features - 
they are perfectly continuous, finite and station- 
ary. Their shape depends on the assumptions we 
make about disk vertical structure. It is possible 
that for more realistic models they would be even 
more profound. 
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Fig. 5. — The non-monotonic sections of the ra- 
dial velocity and surface density profiles in an ac- 
cretion disk surrounding a 9.4 M© BH with a = 
0.1. The solid lines are for mass accretion rates 
between 0.1 and 0.9rriEdd- 

The local minima in the radial velocity profiles 
occur in some particular range of mass accretion 
rates. The minimal accretion rates for which the 
non-monotonic features occur correspond to the 
appearance of radiation pressure dominated disk 
regions. The maximal values are determined by 
the rate of advection which significantly changes 
the disk structure. Fig. [6] presents relation be- 
tween the mass accretion rates for which we ob- 
serve non-monotonic sections in the radial veloc- 
ity profiles and the BH spin. For non-rotating BH 
the local minima appear for accretion rates be- 
tween 0.1 and 0.7rriEdd while for a highly spin- 
ning BH (a* = 0.9) for accretion rates range 
0.025 4- 0.2rriEdd- The solid line presents mass ac- 
cretion rates for which the minimum is the deepest 
(compare Fig. [5]). In general the minima are most 
profound for mass accretion rates lying almost per- 
fectly in between the limiting values i.e. 0.38 for 
a* = and 0.11 for a* = 0.9. These deepest min- 
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ima occur at radii plotted in the bottom panel of 
Fig. [6l They lie just outside the innermost stable 
orbit for a given BH spin. 



■6 0.6 




Fig. 6. — The shaded area at the upper panel 
presents mass accretions rates for which the non- 
monotonic sections of the radial velocity profile 
(see t l4.2p appear at a given BH spin. These fea- 
tures are the most profound at the mass accretion 
rate given by the solid line. The bottom panel 
presents location of the most profound local mini- 
mum in the radial velocity profile as a function of 
BH spin. Location of the innermost stable orbit is 
also presented. 

The bottom left panel of Fig. |4] presents the 
H/r ratio for a few values of the accretion rate. 
For disks with the lowest rh, which are gas pressure 
dominated only, the H/r ratio is almost constant 
down to ~ 10M and then rapidly decreases. For 
higher accretion rates there is a radiation pressure 
dominated inner region of a disk which results in 
an increase of H/r ratio. 

The central temperature profile dependence on 
the mass accretion rate is presented in the bottom 
right panel of Fig. |H The exponents of the tem- 



perature profiles outside ~ 10M depend on the 
gas to radiation pressure ratio. The temperature 
in the gas dominated regions increase more rapidly 
with decreasing radius than in the inner radiation 
pressure dominated part of a disk. As in the case 
of the flux profiles (Fig. [3]) the location of the tem- 
perature maximum moves closer to the ISCO with 
increasing mass accretion rate. In the plunging re- 
gion the temperature drops rapidly and the disk 
becomes gas pressure dominated. 

4.3. Angular momentum profiles 

In Fig. [7] we present profiles of the specific an- 
gular momentum I — u^/ut of the accreted gas 
for few different mass accretion rates in case of 
a non- and highly spinning BHs. For low and 
moderate mass accretion rates the specific angu- 
lar momentum profiles are very close to the Ke- 
plerian distribution. They differ significantly for 
high mass accretion rates. As one could expect 
for an accretion disk ( Abramowicz fc Straubl feoOS) 
the flow is sub-Keplerian at large radii (the higher 
mass accretion rate the more significant deviation 
from the Keplerian flow) . There is a transition ra- 
dius (r center) where the specific angular momen- 
tum crosses the Keplerian profile. This radius cor- 
responds to the location of the pressure maximum. 
Inside r cen ter the angular momentum is superke- 
plerian and crosses again the Keplerian profile at 
ri n interpreted as the disk inner edge. The spe- 
cific angular momentum profiles approach the val- 
ues corresponding to the angular momentum at 
the inner edge of a disk (£j n ). 
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Fig. 7. — Accretion disk specific angular mo- 
mentum profiles for non-spinning (a* = 0.0) and 
rapidly rotating (a* = 0.9) BHs. For each value 
profiles for two mass accretion rates are presented: 
m = 0.1 and rh = B.OriiEdd- 

In Fig. [8] we present the dependence of the 
location of disk characteristic radii (rj n , r$ and 
reenter) on the mass accretion rate for two values 
of BH angular momentum. For low mass accretion 
rates (rh < O.&rriEdd for a* — and rh < 0.2rhEdd 
for a* — 0.9) the location of the disk inner edge 
(ri n ) coincides with the location of the sonic point 
(rs) and is almost independent of the accretion 
rate. In this regime the disk inner edge is located 
very close to the ISCO (marked with dotted lines). 
For higher mass accretion rates both rj n and rs 
move closer to the horizon. Initially, the sonic 
point moves inward slower and it does not coincide 
with the inner edge of the disk. For the highest 
mass accretion rates they are again located close 
to each other. The location of the pressure maxi- 
mum (r C enter) is roughly independent of the accre- 
tion rate for low accretion rates. For rh > 0.6iriEdd 
in case of a non-rotating BH and rh > 0.2mEdd for 
a* = 0.9 it moves slightly inward reaching ~ 10M 
for a non-rotating and ~ AM for a highly-spinning 
BH at rh = lOrriEdd- 

A more detailed discussion of several astrophys- 
ically interesting issues connected to the location 
of the inner edge of slim disks, that follow from 
the transonic models calculated here, will be soon 
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Fig. 8. — Locations of the characteristic points 
of a disk for non-spinning (a* = 0.0, thick lines) 
and rapidly rotating (a* = 0.9, thin lines) BHs for 
a = 0.1. Dashed lines mark location of the inner 
edge of the disk, solid lines stand for the sonic 
point while dot-dashed lines are for the location 
of the pressure maximum. Dotted lines denote 
location of the innermost stable circular orbits. 



4.4. Radiative efficiency 

In Fig. [9] we plot the radiative efficiency of ac- 
cretion, defined as: 



r) 



L 1 L/L E dd 
mc 2 16 rn/rhEdd 



(13) 



versus mass accretion rate for a few values of the 
spin. As discussed by e .g. iKozlowski et al 



BH 



(197 



or 



Jaroszyhski et al.l (198 1 



for accretion 
disks that are thin on their inner edge, the ef- 
ficiency of accretion rj can be approximated by 
r\ = l — ut(ri n ). As it was discussed in the previous 
paragraph, when the accretion rate is very small, 
the inner edge of the disk coincides with the loca- 
tion of the ISCO (Fig. \8§. Thus, for very small 
accretion rates, the efficiency is constant (does not 
depend on the accretion rate) , and the total lumi- 
nosity is proportional to the accretion rate. 

It was found by Jaroszyhski et al.l (1980). 



Abramowicz et al. ( 198g ). |Paczyhskil l|l982h jpaczyhskil 
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(199^) and several other authors, that when accre- 
tion rate increases, the inner edge of a disk moves 
from the location of ISCO to the location of the 
marginally bound orbit r mb (cf. Fig. [8]). Because 
Ut(i~mb) — 1, the efficiency at tmb is formally zero, 
and therefore, as found in the above quoted pa- 
pers, when accretion rate increases, luminosity is 
not proportional to the accretion rate, but grows 
more slowly. We illustrate the drop in efficiency 
in Fig. EE 

A more detailed discussion of the efficiency of 
accretion at h igh accretion rates may be found in 
the review by lAbramowicz ( 2005 ). 




accretion rate [rn Edd ] 

Fig. 9. — The top panel present the total lumi- 
nosity of an accretion disk for different values of 
BH spin as a function of mass accretion rate. The 
bottom panel presents the efficiency parameter r) 
(see 30). 



4.5. Alpha dependence 

The a formalism I Shakura fe Sunvaev 1973h 



gives only the upper limit for the value of a pa- 
rameter which has to fullfil the condition a < 1. 
This restriction comes from the assumption that 
turbulent elements of size smaller that disk thick- 
ness and moving at velocities smaller than speed 
of sound waves are the source of viscosity. The 
usual approach is to set a to some particular con- 
stant value. In this paper we choose a = 0.1 as 
our standard assumption. 



Fortunately, as Shakura fc Sunvaev 1 19731 ) 
have proven, the outcoming flux of energy does 
not depend on a in the case of radiatively efficient 
accretion disks. Whether this is true or not for 
disks with advection has to be checked numeri- 
cally. In Fig. [10] we present our solutions for a few 
different values of the a parameter and for two 
mass accretion rates. The surface density pro- 
files depend strongly on a - to generate the same 
amount of energy a model with lower value of a 
requires higher value of the vertically integrated 
pressure which corresponds to higher column den- 
sity. According to the bottom panels of Fig. [TO] 
the flux profiles are insensitive to different values 
of a at all radii outside the ISCO even for high 
mass accretion rates. The differences in flux pro- 
files inside the ISCO do not influence the emergent 
spectra as the emission in this region is negligible. 
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Fig. 10. — Surface density and flux profiles de- 
pendence on the value of the a parameter for two 
mass accretion rates m = O.lOmEdd (left panels) 
and l.OOmEdd (right panels). BH mass is 9.4 M©. 

In Fig. [TTJ we present the location of the sonic 
point for different values of the a parameter and 
different mass accretion rates for a non-rotating 
BH. For very low accretion rates the sonic point 
is located close to the ISCO independently of the 
value of the visosity parameter. For high mass ac- 
cretion rates and low viscosities the sonic point, 
as the theory of thick accretion disk predicts (e.g. 
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Jaroszvriski et al. 1980l : Paczvriskil 19821 ). moves 
inward towards the marginally bound orbit. For 
the highest values of the a parameter (a > 0.2) 
the behaviour is opposite - the sonic point moves 
outside the ISCO with increasing mass accretio n 
rates (compare Fig. 11 in Abramowicz et alI T988). 
Such solutions are interpreted as Bondi -like and 
are not unique I Muchotrzeb-Czernylll986f ). There- 
fore the dotted lines in Fig. QJ] should be inter- 
preted as profiles describing approximate mean lo- 
cations of the sonic point in this regime. 




Fig. 11. — Location of the sonic point as a func- 
tion of mass accretion rate for different values of 
the a parameter: a = 0.001, 0.01, 0.1, 0.2, 0.3, 
and 0.4 for a non-rotating BH. Dotted lines (for 
highest values of a) denote approximate mean lo- 
cation of the sonic point as the solutions in this 
regime are not unique (see t )4.5p 

5. Parameter study 

The main goal of this work is to provide a spec- 
trum of slim disk solutions in a wide range of pa- 
rameters applicable to stellar mass X-ray binaries. 
The parameter space is spanned by four values 
defining the system properites: BH mass, accre- 
tion rate, BH angular momentum and the value of 
a. In the standard Novikov & Thorne approach 
the flux profile (most important for spectral con- 
tinuum fitting) does not depend on the value of a 
and is inversely proportional (for given r/M and 



in/iTiEdd) to the BH mass. Our results show that 
taking advection into account does not change 
these relations, e.g. the flux profile hardly depends 
on a ( §4.5p neither the advection rate (cf. Fig. [2]) 
does depend on the BH mass. Therefore, the ma- 
jor properties of an accretion disk depend mostly 
on two parameters: the mass accretion rate rh and 
the dimensionless spin parameter a* = a/M. In 
this section we point out the existence and sum- 
marize the properties of three regimes in rh corre- 
sponding to three different bra nches of the M(S) 
curve described in details in lAbramowicz et ; 
( 19881 ). The mass accretion rates limiting them 



depend on BH spin and coincide with accretion 
rate limits for the non-monotonic features in the 
radial velocity profiles presented in the top pannel 
of Fig. El 

(i) low accretion rates — Disk is gas dominated 
only. The radiation pressure supported inner re- 
gion does not occur. The radial velocity profile is 
monotonic. Disk is radiatively efficient {rj « 1 for 
a* = and rj I for a* — > 1). The advection 
is not significant. The flux profile agrees with the 
Novikov & Thorne solution. The angular momen- 
tum profile is almost Keplerian. Both r in and rs 
are very close to the ISCO. 

(ii) moderate accretion rates — Radiation pres- 
sure dominates the inner region of a disk. The 
radial velocity and surface density profiles are no 
longer monotonic. Disk is not radiatively efficient: 
r) drops down with increasing accretion rate. The 
advection becomes more significant what results in 
a shift of the flux profile maximum with respect to 
the radiatively efficient case. Rotation is slightly 
sub-Keplerian outside r cen t e r and super-Keplerian 
between r center and ri n , ri n and rs almost coin- 
cide and are located close to the ISCO. 

(iii) high accretion rates — The radiation pres- 
sure dominated region extends upto large radii. 
The advection becomes dominant significantly 
shifting the flux maximum inward. The radial 
velocity profile is monotonic. The angular mo- 
mentum profile significantly deviates from the 
Keplerian distribution. Disk is radiatively inef- 
ficient (rj — » 0). ri n and rs are located inside the 
ISCO. 
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6. Prospects for future work 

Models of slim accretion disks calculated here 
are needed in several astrophysical applications 
that are at present being worked by us in a frame- 
work of a larger research project (connected to my 
Ph.D. work) that includes: 

(1) Combining the transonic radial descrip- 
tion of the slim disk structure used here with 
a more accurate treatment of the vertical struc- 
ture by using ideas and methods developed to 
solve the vertical radiation transfer in thin accre- 



tion d i sks by I Hameu rv et alJ (I1998T) : iDavis et al . 
l)2005T) : lR6zanska fc Madejl (|2008h and others. We 
are working on a general numerical code that 
would self-consistently combine radial and vertical 
structure calculations. The major issue here is to 
calculate the location of the effective photosphere 
of a slim disk. We have already tested a simplified 
version of the code and examined a thin, radia- 
tion pressure supported disk, confirming that in 
this case location of the effective photosphere only 
weakly depends on the (unknown and therefore ad 



hoc a ssumed) vertical dissipation IjSadowski et al 
l2009h . 

(2) Calculating the "observed" X-ray continuum 
spectra of slim disks by ray-tracing photon trajec- 
tories (in the Kerr geometry) all the way from an 
emission place at the disk effective photosphere to 
a distant observer. We use a par ticular versio n of 
the ray-tracing code developed by Bursa ( 2006T) . It 
fully includes all special and general relativistic ef- 
fects. An extensive catalog of t he calculated spec- 
tra w ill be published elsewhere ([Bursa &: Sqdowski 
2009) and used in collaboration with the Harvard 
group l Abramowicz et al. 20091 ) to improve the 
black hole spi n estimates by the spectral X-ray fit- 
ting ( as in e.g. Shafee et al. 20061 : Middleton et al 



2006: iGou et al.l l2009l. Improvements follow fro 



including effects of advection (they are relevant for 
higher accretion rates) and from a more accurate 
calculation of the location of the effective photo- 
sphere (relevant at higher inclinations). 

(3) The stationary slim disk models calculated 
here provide the initial conditions needed in nu- 
merical simulations of non-stationary slim disks 
in Kerr geometry that we perform in collabora- 
tion with the Xiamen University group. The inter- 
esting astrophysical issue here is a possible limit- 
cycle behavior. Also in this case we are calculating 



the observed appearance of the disk using the ray- 
tracing from the effective photosphere. 

7. Conclusions 

In this paper we presented a numerical method 
used to solve equations describing slim accretion 
disks. The most important assumptions we made 
were as following: we assumed stationarity and 
axis symmetry, we neglected the angular momen- 
tum flux taken away by radiation and the radial 
flux of radiation, we assumed large optical depth, 
used the vertical equilibr ium formula derived by 
Abramowicz et alJ (| 1997T I and we allowed for ad- 
vective flux of energy. Slim disk equations were 
reduced to a two dimensional two boundary value 
problem. W e solved it applying relaxation method 
(<Pressll2002T ). The solutions are presented in de- 
tails in |H For moderate and high mass accre- 
tion rates we observe large amount of advection 
which significantly changes the emergent flux pro- 
files. For some particular range of accretion rates 
non-monotonic features in disk structure appear. 
We conclude they could be even more profound 
under different assumptions about disk vertical 
structure. The applications of the relativistic slim 
disk solutions presented here are discussed in $6l 
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A. The manual for tabulated slim disk solutions available online at 
http : //users . camk . edu . pi/ as/slimdisk 

The solutions of the relativistic slim accretion disk model presented and discussed in this paper have been 
available online at http://users.camk.edu.pl/as/slimdisksince October 2009. In this section we present 
a brief description of the format they are provided in as well as for the included C interpolation routines. 

There are two sets of files available: one for 10 M Q BH only (weighting ca. 170 MB) and the other for a 
few BH masses with mass interpolation included (ca. 1 GB). The former is enough for ray-tracing purposes 
where the effective temperature and disk thickness profiles are the only needed while the latter for more 
general applications where other quantities, depending on BH mass in more complicated ways, are involved. 

The disk profiles are provided in the following subfolders: data-dsi/ in case of the single mass set of 
solutions and data-dsi-mN/ in case of the full set of disk models where N stands for different BH masses 
and can take the following values: 5, 10, 20, 30, 50 and 100 M Q . Each folder contains a number of 
soltt.N.M.dat files describing accretion disks for given set of the input parameters. These files are indexed 
in the res.mamdot.dat file containing 8 columns with the following meaning: 

1. BH mass [M Q ], 

2. disk luminosity defined as: 2nrF em dr in units of the Eddington luminosity: Lsdd = 47T ^ Mc = 
1.25 xlO 38 -^ erg-s" 1 , 

3. mass accretion rate in units of the critical mass accretion rate defined in Eq. [TTJ 

4. dimensionless spin parameter a* = GI ^2/ C , 

5. location of the sonic point in units of 2M, 

6. value of the angular momentum at the horizon (the eigen value of the transsonic solution) in units of 

2M, 

7. N file index, 

8. M file index. 

Each soltt .N.M.dat file contains 12 columns in the following order. All values but radius and temperature 
are given in G = c = 1 units (conversion factors are given in Tab. [J). 



1. 


radius [2M], 


2. 


radial velocity (as defined in Sect. [2]) 


3. 


central temperature [K], 


4. 


angular momentum, 


5. 


surface density, 


6. 


H/r ratio, 


7. 


vertically integrated pressure, 


8. 


radiation to gas pressure ratio, 


9. 


flux emitted (both sides), 


10. 


flux adverted, 



13 



11. specific angular momentum (u^/ut), 

12. undefined. 

The emitted flux F em is related to the disk effective temperature T e // by the following formula: F em = 
2crT^f where a = 1.56 • 1CP 60 m _2 K _1 is the Stefan-Boltzmann constant in G = c = 1 units. 

Table 1: Conversion factors from cgs to geometrical units 



velocity 


1 cm/s = 3.34 • 1CT 11 


angular momentum 


1 g cm 2 /s = 2.48 • 1CT 43 m 2 


surface density 


1 g/cm 2 = 7.42- lCT 27 1/m 


vertically integrated pressure 


1 Ba cm = 8.26 • 1CT 48 1/m 


flux 


1 crg/cm 2 /s = 2.76 • 10~ bb 1/m 2 



Disk solutions can be accessed directly or through provided interpolation routines. Two C files are 
attached to the archives: dsi.c and draw_dsi.c. The former contains the interpolation routines and is 
designed to be easily included in any C code requiring instant access to disk solutions at any combination 
of the input parameters (mass accretion rate or luminosity, BH spin, radius and BH mass). The latter 
extracts given disk profiles for chosen input parameters, prints them to files and draw the solutions to 
.GIF and PostScript files using GnuPlot. Compilation of these routines requires GNU Scientific Library 
(http://www.gnu.org/software/gsl). The header file dsi.h as well as an exemplary Makefile are also 
included. 

The interpolation routines take the same input parameters for both archives. The only difference is that 
the routines included in the smaller archive neglect the mass and always returns solutions for 10 Mq. 
To include slim disk interpolation into a C code one should do the following: 

1. include dsi .h header file, 

2. set up the interpolation by calling dsi_setup(MdotLum, a, M, alpha, MdotLumFlag, fermin, 
fermax) with the following arguments: 

• MdotLum - mass accretion rate or disk luminosity in Eddington units depending on MdotLumFlag, 

• a - dimensionless spin parameter, 

• M - BH mass [ Mq], 

• alpha - neglected, 

• MdotLumFlag - for disk luminosity as the 1st argument, 1 for mass accretion rate, 

• rmin - minimal radius of the tabulated solution is returned, 

• rmax - maximal radius of the tabulated solution. 

The routine returns if the input parameters lie inside the grid of tabulated solutions and 1 if go 
outside the grid at least in one dimension. 

3. call dsi_eval(r, N) to evaluate a given parameter at radius r (in units of 2M). The meaning of N is 
the following: 
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absolute value of the radial velocity 

1 central temperature 

2 specific angular momentum 

3 surface density 

4 H/r ratio 

5 total vertically integrated pressure 

6 radiation to gas pressure ratio 

7 flux emitted 

8 flux advected 

9 angular momentum 

All values are returned in geometrical units as discussed above. 
4. set off interpolation for given disk parameters by calling dsi_setof f () . 

The draw_dsi.c file contains additional stand alone code which uses the routines provided by dsi.c to 
print out and draw disk solutions for given system parameters. The code takes input in the following format: 
./draw_dsi MBH N ml If a 7„f [a 7„f ml '/.f] [a/ml °/ D f] filename 
where the meaning of the arguments is following: 

• MBH - BH mass [M ], 

• N - index of the required quantity (as given above) , 

• ml - disk luminosity (if negative) or mass accretion rate (if positive) , 

• a - dimensionless BH spin, 

• filename - name of output file for .gif and .ps figures. 
Examples: 

./draw_dsi.out 10 3 ml .3a0a .6a .9 figl 

- plots surface density for M=10M Q and the following (mdot, a) sets: (.3, 0), (.3, .6), (.3, .9) 
./draw_dsi.out 10 7 a ml - . 1 ml -.2 ml -.3 figl 

- plots flux for M=10M Q and the following (luminosity, a) sets: (.1, 0), (.2, 0), (.3, 0) 
./draw_dsi.out 10 7 a ml .6a .9 ml .3a .99 ml .1 figl 

- plots flux for M=10M Q and the following (mdot, a) sets: (.6, 0), (.3, .9), (.1, .99) 

Profiles of the required quantity for each set of disk parameters are printed out to out_n.dat files each 
containing two columns with radius in the first one and profile of the required value in the second column. 

The slim disk solutions are tabulated for a = 0.1, < a* < 0.99 and M < 500MEdd- 
The author will appreciate any comments and bug reports sent to as@camk.edu.pl. 



15 



